# rank the perturbations

args = commandArgs(TRUE)
cutoff = args[1]

infile = "cutoff%s-perturbation_%s_interfaceNode.txt"

done4 = FALSE
if(file.exists(sprintf(infile,cutoff,4))) done4 = TRUE

p1 = read.table(sprintf(infile,cutoff,1),header=TRUE,row.names=1, stringsAsFactors=FALSE,sep="\t")
p2 = read.table(sprintf(infile,cutoff,2),header=TRUE,row.names=1,stringsAsFactors=FALSE,sep="\t")
p3 = read.table(sprintf(infile,cutoff,3),header=TRUE,row.names=1,stringsAsFactors=FALSE,sep="\t")
if (done4) p4 = read.table(sprintf(infile,cutoff,4),header=TRUE,row.names=1,stringsAsFactors=FALSE,sep="\t")

unp = p1["Unperturbed",]

p1=p1[!(rownames(p1) == "Unperturbed"),]
p2=p2[!(rownames(p2) == "Unperturbed"),]
p3=p3[!(rownames(p3) == "Unperturbed"),]
if (done4) p4=p4[!(rownames(p4) == "Unperturbed"),]

perts = rbind(p1,p2,p3)
if (done4) perts = rbind(perts,p4)
perts = perts [perts$Simulation_count < 1001,] # remove perturbations that didn't converge in any attractor

# number of perturbed genes == number of "_"
library(stringr)
perts$Perturbed_nodes = str_count(rownames(perts), "_")
perts$avg_flipped_genes = perts$Flipped_genes / perts$Perturbed_nodes

perts = perts[order(perts$avg_flipped_genes, decreasing = TRUE),]

write.table(perts, sprintf("cutoff%s_all_ranked_perturbations.txt",cutoff), sep="\t")

best = 100-min(perts$Match_fraction_initial)

print(sprintf("Perturbed combinations of 1-%s TFs",ifelse(done4,4,3)))
print(sprintf("Max flipping is %s%%",best))
if (best < 40) print("NOT FLIPPING ENOUGH")